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We present a detailed analysis, based on the Forward Flux Sampling (FFS) sim- 
ulation method, of the switching dynamics and stability of two models of genetic 
toggle switches, consisting of two mutually-repressing genes encoding transcription 
factors (TFs); in one model (the exclusive switch), the two transcription factors mu- 
tually exclude each other's binding, while in the other model (general switch) the 
two transcription factors can bind simultaneously to the shared operator region. We 
assess the role of two pairs of reactions that influence the stability of these switches: 
TF-TF homodimerisation and TF-DNA association/dissociation. In both cases, the 
switch flipping rate increases with the rate of TF dimerisation, while it decreases 
with the rate of TF-operator binding. We factorise the flipping rate k into the prod- 
uct of the probability p{q*) of finding the system at the dividing surface (separatrix) 
between the two stable states, and a kinetic prefactor R. In the case of the exclusive 
switch, the rate of TF-operator binding affects both p{q*) and i?, while the rate of TF 
dimerisation affects only R. The general switch displays a higher flipping rate than 
the exclusive switch, and both TF-operator binding and TF dimerisation affect k, R 
and p{q*)- To elucidate this, we analyse the transition state ensemble (TSE). For the 
exclusive switch, the TSE is strongly affected by the rate of TF-operator binding, 
but unaffected by varying the rate of TF-TF binding. Thus, varying the rate of TF- 
operator binding can drastically change the pathway of switching, while changing 
the rate of dimerisation changes the switching rate without altering the mechanism. 
The switching pathways of the general switch are highly robust to changes in the rate 
constants of both TF-operator and TF-TF binding, even though these rate constants 
do affect the flipping rate; this feature is unique for non-equilibrium systems. 

PACS numbers: 



I. INTRODUCTION 

Biochemical networks with multiple stable states are omnipresent in living cells. Multi- 
stability can provide cellular memory, it can enhance the sharpness of the response to intra- 
and extracellular signals, it can make the cell robust against biochemical noise, and it allows 
cells to differentiate into distinct cell types. Although a multistable biochemical network 
can flip between alternative states due to random fluctuations ("noise"), in many cases the 
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states are very stable and the network typically only switches from one state to the next 
under the influence of an external signal A key question, therefore, in understanding 
multistable biochemical networks is what controls the stability of the steady states. To 
answer this question, we have to elucidate the pathways of switching between steady states. 
Switching events are, however, intrinsically difficult to study experimentally, because the 
switching event itself can be much faster than the typical life time of the steady state. 
Computer simulations are a valuable tool for studying biochemical networks, especially for 
rare processes such as switching. However, precisely because such events are rare, special 
techniques are required to simulate them. One such technique is Forward Flux Sampling 
(FFS) 0,1,3, and in this paper, we use FFS in combination with committor distributions 
to analyse in detail the effect of two important sources of fluctuations — transcription fac- 
tor dimerisation and transcription factor-DNA binding — on the flipping rate and switching 
pathways of two models of bistable genetic toggle switches. We hope that this analysis may 
serve as a paradigm for studying multistable biochemical networks as well as other rare 
events in non-equilibrium systems. 

If a biochemical network is bistable, with two stable states A and B, respectively, then 
it will show a bimodal steady-state probability distribution, p{q), of some order parameter 
q. This order parameter can be the concentration of a species, or a combination of the 
concentrations of a number of species. It is usually interpreted as a reaction coordinate 
that measures the progress of the 'reaction' from state A to B. Recently, such bimodal 
distributions have been measured experimentally for biochemical networks jl, @, 0]- These 
distributions are potentially useful, because they are linked to the rate of switching from one 
state to the other. In particular, we have recently shown [sl that not only for equilibrium 
systems, but also for systems that are out of equilibrium such as biochemical networks, the 
rate of switching from state A to state B, kAB, can be written as the product of two factors: 

kAB = Rpiq*)- (1) 

Here, q* denotes the location of the dividing surface, the separatrix jl, 0], which separates 
the two states A and B. The above relation is useful because it shows that the rate of 
switching from one steady state to the next, is given by the probability p{q*) of being at the 
dividing surface times a kinetic prefactor R that describes the average flux of trajectories 
crossing the dividing surface. However, while the rate constant kAB does not depend on 
the choice of the order parameter q as long as it connects states A and B, both p{q*) and 
R do depend on the choice of q. If q is the "true" reaction co-ordinate that accurately 
describes the switching process, then q* corresponds to the transition state and p{q*) and R 
provide accurate measures for the probability of being at the transition state and the flux 
of trajectories leaving the transition state for state B [13]. A key issue in analysing rare 
events in general is therefore identifying the reaction coordinate q that accurately describes 
the progress of the transition. 

FFS can be used to compute kAB, p{q*), R, and to generate members of the transition 
path ensemble 0, 0, 0, [HI . To identify the reaction coordinate, the transition paths can be 
analysed using committor distributions; this approach is commonly applied in the field of 
soft-condensed matter physics [lo|, and we have recently demonstrated how it can be used 
to analyse the transition pathways of biochemical switches j2[ . Each configuration x of our 
system has a commitment probability or "committor" , Pe (x) ■ This is the probability that a 
trajectory, fired at random from that configuration, will reach state B before state A. Given 
Pb{x), we can define the "transition state ensemble" (TSE) [13], which is the collection 
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of configurations along tlie reaction pathways which have committor value Pb{x) = 0.5. 
We can extract TSE configurations from our switching pathways by computing committor 
values for all the points along the pathways and selecting those points (several per path) 
with Pb = 0.5. We then try to find an order parameter (or combination of order parameters) 
that accurately describes these TSE configurations. To test likely order parameters, one can 
compute the probability distribution for the order parameter for the TSE configurations 
[2I, Poorly chosen order parameters will show a broad or even bimodal distribution 13|. 
while good order parameters will show a narrow distribution of values in the TSE. 

In this paper, we apply FFS to study two different models of genetic toggle switches, 



consisting of two genes A and B that mutually repress each other [2, la, [IJ, ll5|, [ig, ll7|, ll8|, 
19 . 20| . The genes encode transcription factors (TP) A and B respectively. These can form 
homodimers, in which form they can bind to a regulatory region of the DNA (represented by 
an operator site O) and regulate transcription. The dimer A2 represses the transcription of 
B when bound to O and vice versa (see Pig. [T]). In the first model, called the exclusive switch 
[sl, [13] the dimers of the two species are not allowed to simultaneously bind to the operator; in 
the second model, called the general switch, the operator can bind both types of homodimers 



at the same time [8|, \17j\. Both switches have one stable state in which A is abundant, and 
B scarce, and another stable state in which B is abundant and A scarce. We simulate the 
switch using the Gillespie algorithm [2]?], in combination with PPS. The Gillespie Algorithm 
is a widely used and efficient Kinetic Monte Carlo scheme 22[] for chemical reactions, which 
generates trajectories in correspondence with the chemical master equation. 

Switching events are driven by random fluctuations. Key fluctuation sources in this net- 
work are TP-TP and TP-DNA association and dissociation reactions. By varying the rates of 
these reactions, while keeping their equilibrium constants fixed, we can vary independently 
the time scales and hence the correlation times of these fluctuations. The correlation times 
of fluctuations are important, since they determine the extent to which the fluctuations 



propagate in the network [23|, |2J, |25 



We therefore begin by calculating how the stability of both switches varies with the rate 
of TP-TP association and dissociation (dimerisation), and with the rate of TP-operator 
association and dissociation (operator binding). We vary the association and dissociation 
rates together, keeping their ratio {i.e. the equilibrium constant) unchanged. The switching 
rate is strongly affected: for both models, decreases as the rate of operator binding 
increases, and increases as the rate of dimerisation increases. Analysing the effects on p{q*) 
and R, we find that the two models behave differently: for the exclusive switch, the rate of 
operator binding changes both p{q*) and R, while the rate of dimerisation affects only R; 
for the general switch, the rate of dimerisation affects both p{q*) and R, while the rate of 
operator binding predominantly changes p{q*). 

We then show that the effect of TP-TP and TP-DNA fluctuations on k, R, and p{q*) 
can be understood by elucidating the switching mechanism using commmittor distributions. 
We find that for the exclusive switch the difference in total copy number number of the 
two species is not a complete reaction coordinate: the state of the operator is also an 
important factor in determining the committor value |^]. In contrast, we find little evidence 
that dimerisation is an important ingredient of the reaction coordinate. This explains why 
the rate of operator binding affects both the probability of being at the separatrix and 
the kinetic prefactor, while dimerisation only affects the kinetic prefactor. Por the general 
switch, the situation is markedly different: the switching mechanism is highly robust to 
changes in both the rate of operator binding and the rate of dimerisation. Hence, changing 
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these rate constants does not change the route the switching pathways take in state space, 
yet does affect the ffipping rate. This is a manifestation of the fact that this is a non- 
equihbrium system — in an equihbrium system the switching rate cannot be changed without 
changing the switching pathways. The general imphcation of this observation is that in 
order to understand the stabihty of biochemical switches, we need to understand not only 
the composition of the transition state ensemble, but also the dynamics of the transition 
paths. 

In the next Section, we describe the model genetic switches in more detail. In the sub- 
sequent Section, we briefly discuss the FFS technique. We then present the results on the 
switching rate, the kinetic prefactor and the probability of being at the separatrix for both 
models, showing how they depend on the rates of operator binding and of dimerisation. In 
the next Sections, we discuss switching pathways and reaction coordinates first for the ex- 
clusive switch, and then for the general switch. We end with a discussion of the implications 
of our findings for the modelling of multistable biochemical networks and the study of rare 
events in other non-equilibrium systems. 



II. MODELS: THE EXCLUSIVE AND THE GENERAL SWITCH 



We consider a genetic toggle switch consisting of two genes, each of which represses 
the other 0, [s], [ij, 17, [Tsl, 26]. A switch of this kind has been constructed and shown 
to be bistable in E. coli l5[. We study both the general switch and the exclusive switch, 
introduced by one of us [8|, [iTJ. The general switch is represented by the following set of 
reactions [l, S, 17]: 



A + A^^ A2, 

fcb 


fcf 

B + B 4 B2 

feb 


(2a) 


+ A2 OA,, 

feoff 


+ B2 OB2 

feoff 


(2b) 


0A2 + B2 OA2B2, 

feofl 


OB2 + A2 OA2B2 

feoff 


(2c) 




O'^-^'O + B 


(2d) 


OAa'^'OAa + A, 


OB2'^^'OB2 + B 


(2e) 


AA0, 


B^0 


(2f) 



In this reaction scheme, O represents a DNA regulatory sequence adjacent to two divergently 
transcribed genes A and B. These code respectively for proteins A and B, as shown in Fig. 
[U Genes A and B can each randomly produce proteins with the same rate, but whether they 
do so depends on the state of the operator O. Proteins A and B can each form a homodimer 
that can bind to the operator. When an A2 dimer is bound to O, the production of B 
is blocked, and likewise, when a B2 dimer is bound to O, the production of A is blocked. 
When both dimers are bound to the operator, no protein can be produced. Proteins can also 
vanish (in the monomer form), modelling degradation and dilution in a cell. This model can 
be modified by removing reactions (^): in this case, transcription factors mutually exclude 
each other's binding to the operator. We refer to the switch described by the whole set of 
chemical reactions ([2]) as the "general switch"; the "exclusive switch" consists of the same 
set of reactions, except for reactions 1^). 
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We have assumed in this model that transcription, translation and protein folding can 
be modelled as single Poisson processes, neglecting the many substeps that lead to the 
production of a protein. Ref. [8l] discusses the effects of both shot noise and fluctuations 
in the number of proteins produced per mRNA transcript on the switch stability. We also 
note here that while mean-field analysis predicts that cooperative binding of the TPs to the 
DNA is essential for bistability [2^, it has recently been demonstrated that bistability can 
be achieved without cooperative binding when the discrete nature of the reactants is taken 
into account {U 



We choose as the unit of time for our simulations, and we use the volume of the 
system, V, as the unit of volume. In this paper, we will use the following "baseline" set 
of parameters: k{ = 5kprodV, kh = Sfcprod (so that = fcb/fcf = 1/^), ^on = Sfcprod^, 
kos = /Cprod (so that = kos/kon = 1/{5V)), fi = 0.3/cprod- These parameters are chosen 
to be representative of typical cellular values, as discussed in Section |Vl Por simplicity, the 
model switches are completely symmetrical - rate constants for equivalent reactions involving 
A and B are the same. The mean field analysis performed in demonstrates analytically for 
both systems the existence of three fixed points for the parameter values listed above: two 
symmetrical stable states, one rich in A and the other rich in B, separated by one unstable 
state where the total number of A equals the total number of B. This implies that the system 
can be considered a truly bistable switch. However, while this analysis indicates the regions 
in parameter space where the system is bistable, it cannot predict the switch stability nor 
elucidate the switching pathways. Por this reason, we carry out stochastic Kinetic Monte 
Carlo simulations using the Gillespie algorithm 21, 2^- In previous work, we have shown 



that the switch stability depends strongly on the mean copy number of species A and B [17 
which is given by the ratio of the protein production and decay rates, fcprod/ In this paper, 
we investigate its dependence on the other parameters fcf, k^, kon and kos, which govern key 
sources of fluctuations in the network - TP-TP and TP-DNA association and dissociation 
reactions. 



III. METHOD: FORWARD FLUX SAMPLING 

Conventional simulation methods are ineffective for studying rare events such as the 
flipping of biochemical switches, because the vast majority of the computational effort is 
spent in simulating the uninteresting waiting times in between the events. Por this reason, 
specialised methods are required, and we have recently developed the Porward Plux Sampling 
(PPS) technique 0,1, ||. PPS is well suited to simulating biochemical networks, since unlike 
most other rare event methods [13], it can be used for out-of-equilibrium systems. In this 
paper, we use PPS to calculate rate constants, transition paths and stationary probability 
distribution functions for the model genetic switch. 

To obtain the rate constant kAB for a rare transition between two states A and B, PPS 
exploits the fact that (in steady state) kAB can be written as the product of two factors: 

kAB = ^aPab (3) 

Here, ^a is the number of trajectories that leave state A per unit time, while Pab is the 
conditional probability that these trajectories subsequently reach state B without returning 
to A. An order parameter A must be chosen, which defines states A and B: if A < Aq the 
system is in state A, while it is in state S if A > A^. The parameter A is then used to further 
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define a series of nonintersecting interfaces {Ai, . . . , A„_i}, with Aj < Aj+i, such that any 
trajectory from A to B has to cross all the interfaces {Aq, . . . , A„}, without reaching Aj+i 
before it has crossed Aj. The conditional probability Pab can be written as 

n—l 

PAB = l[P{\^+l\X^), (4) 

i=0 

where P(Aj+i|Aj) is the conditional probability that a trajectory that comes from A and has 
crossed Aj for the first time, will subsequently reach Aj+i before returning to A. Several dif- 
ferent algorithms can be used to calculate $4, P(Ai+i|Ai) and to obtain transition paths; in 
this paper, we have used the original scheme [3, 0] • Briefly, one first performs a conventional 
("brute-force") simulation to compute ^a, which is the number of times that the trajectory 
crosses Aq, coming from A < Aq, divided by the total simulation time. When one of these 
crossings occurs, the configuration of the system is stored, so that this simulation run gener- 
ates not only an estimate for <l>^, but also a collection of points at interface Aq. In the next 
stage, one chooses a point at random from this collection, and fires off a new trajectory from 
this point, which is continued until the system either reaches the next interface Ai or returns 
to state A. If Ai is reached, the system configuration at Ai is stored in a new collection. This 
procedure is repeated a number of times until a sufficiently large number of points at the 
next interface has been generated. An estimate for P(Ai|Ao) is obtained from the number of 
trials which reached Ai, divided by the total number of trials fired from Aq. Starting from 
the new collection of points at Ai, one then repeats this whole procedure to drive the system 
to A2, and so on. Eventually, the system reaches state B, upon which the rate constant 
can be calculated from Eqs.([3]) and (jlj). Furthermore, a (correctly weighted) collection of 
trajectories corresponding to the transition ("transition paths") can be obtained by tracing 
back those trial paths that arrive in B to their origin in A. 

We have recently shown [lH that this procedure can be used to generate not only the 
rate constant and transition paths, but also the stationary distribution p{q), as a function 
of a chosen order parameter (or parameters) q. This is achieved by continuously updating 
a histogram in the parameter q during the trial run procedure, as described in Ref. 
To obtain p{q), histograms for the "forward" {A B) and "backward" {B —>■ A) transition 
must be combined. However, since our model switch is symmetric, the two histograms are 
identical in this case. The parameter q does not have to be the same as A, although in this 
paper we have chosen q = X. 

In FFS, a series of interfaces are used to drive the system over a "barrier" , in a ratchet-like 
manner. The efficiency of the method of course depends on the choice of the order parameter 
A, the positioning of the interfaces, number of trials etc [3]. However, it is important to note 
that A does not have to be the true reaction coordinate for the transition. The choice of A 
does not impose any bias on the system dynamics: transition paths are free to follow any 
possible path between state A and B. The choice of A should not affect the computed rate 
constant, transition paths or p{q). Furthermore, the FFS method does not make a Markovian 
assumption about the transition paths, or any assumptions about the distribution of state 
points at the interfaces {Aq, ■ ■ ■ , A„}: each point at interface i lies on a true dynamical path 
which originates in the initial state A. This turns out to be essential for the model genetic 
switch. 

For the FFS calculations presented in this paper, we have chosen as A parameter the 
difference between the total copy numbers of the two proteins: A = ha — ub, with nx = 
Nx + 2A^X2 + 2A^ox2 the total copy number of species X = A or B in the exclusive switch 
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and nx = Nx + 2A^X2 + "^^0X2 + 2Ai'oA2B2 the total copy number of species X = A or B in 
the general switch. 



IV. RESULTS 

Key sources of fluctuations in in this reaction network are TF-TF and TF-DNA asso- 
ciation and dissociation reactions. We can vary the influence of these fluctuations on the 
network dynamics by changing the rate constants for association and dissociation, keeping 
the equilibrium constant (the ratio of association and dissociation rate constants) fixed, so 
that the macroscopic dynamics of the system remain unchanged. When these reactions 
are fast, fluctuations are short-lived on the timescale of the slower protein production and 
degradation reactions, so that the effect of a fluctuation is lost over just a few produc- 
tion/degradation reactions. However, for slow association-dissociation reactions, fluctua- 
tions in, for example, the ratio of monomers to dimers, can persist over the timescale of 
a series of production/decay reactions, and thus have a strong influence on the dynamics 
of the whole network. In what follows, we first discuss the effects of varying the rates of 
operator binding and dimerisation on the switching rate for both genetic switch models, and 
then, to elucidate these effects, we separately discuss the switching pathways for the two 
cases. 



A. Switching rates 

Figs. [2]A and [2]B show the flipping rate kAB for the exclusive switch as a function of the 
dimerisation rate kf and the operator binding rate /Con, respectively (keeping the dissocia- 
tion constants constant). The results for the general switch are shown in Figs. [3]A and[3l3. 
It is clear that for both switches the two sources of fluctuation have very different effects 
on the stability: while kAB increases with the rate of dimerisation (panels A), it decreases 
with the rate of operator binding (panels B). Thus, fluctuations in the TF-DNA associ- 
ation/dissociation reactions tend to destabilise the switch, whereas (counter- intuitively), 
fluctuations in the TF-TF association/dissociation reactions increase the switch stability. 

To understand the origin of these effects, we factorise the flipping rate kAB into the 
product of the probability of finding the system at the dividing surface p{q*) and a kinetic 
prefactor R, as in Eq. ([T]). Fig. H] shows the steady-state probability distribution p(A) of 
finding the system at a particular value of the order parameter A, for different values of the 
dimerisation and operator binding rate: panel A refers to the exclusive switch and panel B 
to the general switch. These distributions were computed using FFS, as described in Section 
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We first note that both distributions exhibit peaks at A ~ ±27, corresponding to the two 
stable states. Secondly, the locations of the two stable states and the shape of the stationary 
distributions are fairly insensitive to both the rate of dimerisation and the rate of operator 
binding. However, around A = the distributions, especially that of the general switch, 
are much more sensitive to changes in the rate constants. Interestingly, the probability of 
finding the system at the value A = is markedly differently for the two models: while for 
the exclusive switch p(A) exhibits a minimum, representing an unstable steady state for the 
system, in the case of the general switch, the probability distribution has a local maximum. 



indicating the presence of a metastable steady state [8|, [iSj. Finally, we note that for an 
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equilibrium system, fluctuations do not influence the stationary distribution function: the 
effect of k{ and k^n on p{q) is a clear characteristic of the non-equilibrium nature of the 
dynamics in this system. 

From the distribution p{q), we compute the probability of being at the minimum of 
the curve, p{q*)- For the exclusive switch, this point corresponds to the dividing surface 
p(A = 0); we show in Figs. EP and[2p how this quantity varies with k{ and fcon, respectively. 
In the case of the general switch, the transition happens through the metastable state at 
A = 0. However, the rate-limiting step for the flipping is to get to the minimum of p(A), 
which is now located at A ~ ±5. Therefore, for the general switch, p{q*) was computed for 
g* = A = 5; it is shown in Figs. [Sp and[3p. By combining p{q*) with Eq. ([1]), we compute 
the kinetic prefactor R, shown in panels E and F of Figs. [2] and [3|, for the exclusive and 
general switch, respectively. 

We observe that for the exclusive switch p{q*) depends upon the operator binding rate 
(Fig. [2P), but not upon the rate of dimerisation (Fig. [2p), while for the general switch 
p{q*) depends upon both rate constants (Figs. [SPjD). In both models, the kinetic prefactor 
R increases with the rate of dimerisation (Figs. [2]Ej3l3), while it decreases with the rate of 
operator binding (albeit much less so in the general switch; Figs. [2]Fj3F). One might expect 
that a change in p{q*) reflects a change in the location of the switching pathways in state 
space. This would suggest that in the general switch, the switching pathways depend upon 
both the rate of dimerisation and the rate of operator binding, while for the exclusive switch 
the switching mechanism does depend upon the rate of operator binding, but not on the 
rate of dimerisation. While the conclusion for the exclusive switch is correct, that for the 
general switch is not, as we discuss in the next two sections. 

B. Switching pathways - Exclusive switch 

To understand the effects of the operator binding and dimerisation fluctuations on the 
switching rate, we would like to determine what the true reaction coordinate is for the 
switching process and whether it involves these fluctuations. To do this, we need to examine 
the transition paths for switching, which are also generated by FFS. We will focus on three 
sets of parameters: (1) the base-line set, with operator binding rate kon = Sfcprod^ and 
dimerisation rate kf = Sfcprod^; (2) a set with slow dimerisation, k{ = O.lfcprod^, and 
fcon = Sfcprod^; (3) a set with fast operator binding, fcon = SOOfcprod^, kf = 5kpj-odV. As 
above, in all cases the dissociation rates are scaled such that the equilibrium constants 
remain constant: = k\,/ki = 1/V and = kos/kon = 1/(5^)- In this section, we 
discuss the exclusive switch, while the next section focusses on the general switch. 

To analyse the progress of the system as it flips from one state to the other, we have 
averaged the switching trajectories in the Pb ensemble. The committor Pb{x) is the prob- 
ability that a trajectory propagated at random from conflguration x reaches state B before 
state A. The Pb ensemble is formed by those conflgurations x that have the same value of 
Pb', {Q{x))pg thus denotes the average of a quantity Q{x) in the ensemble of conflgurations 
X with the same value of Pb- Given an ensemble of switching paths obtained with FFS, we 
can harvest conflgurations x with the same value of Pb- Indeed, each transition path has at 
least one conflguration for every value of Pb- Pb{x) can be used to characterise the progress 
of the transition from A to B — in a sense, it is the true reaction coordinate and our task 
is to identify coordinates that characterise Pb- However, its evaluation is computationally 
very expensive. We have computed Pb for conflgurations in the transition paths that were 
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generated using FFS, by firing 100 test trajectories from eacli configuration. The average 
paths in the Pb ensemble are rather "noisy" , precisely because Pb is a stochastic quantity 
that has to be estimated by a computationally demanding procedure. 

Fig. OA. shows the average switching pathways for the exclusive switch in the nA,n^ 
plane, where ua and ub are the total copy numbers of species A and B, respectively (wa = 
A^A + SA'as + 2A^oA2 for the exclusive switch and ^a = A'a + 2iVA2 + SA'oAa + SiVoAaBa for 
the general switch; similarly for B). Paths are shown for both the A ^ B (solid lines) and 
B ^ A (dashed lines) transitions. Considering the red and black lines in Fig. [5]A., we see 
that the dimerisation rate fcf has little effect on the switching pathways (at least in this 
representation), while on considering the green and black lines, it is clear that the operator 
binding rate fcon does strongly influence the switching pathways, especially in the region of 
the dividing surface, where riA = ?t-b- the pathways shift to lower values of and n-Q when 
fcon is increased. 

Since it appears from Fig. [S]A. that the state of the operator is likely to play an important 
role in the switching mechanism, we plot in Fig. ElA. the probability that the operator 
is bound by a B2 dimer, (A^ob2) as a function of A. Comparing the solid black and red 
lines, we see that changing the rate of dimerisation has indeed little effect on the transition 
paths. In contrast, a comparison of the black and green solid lines shows that changing 
the rate of operator binding has a strong effect on the switching pathways. This indicates 
that operator state fluctuations play an important role in switch flipping ]1, [9| - so that the 
reaction coordinate depends not only on the difference in the number of protein molecules, 
A, but also on which protein is bound to the operator. 

This fact is further illustrated in Fig. [TJ which shows histograms for configurations in 
the TSE of the transition from A to -B; members of the TSE are points along the transition 
paths for which Pb = 0.5. Each panel in Fig. [7] corresponds to a different parameter 
set — the baseline parameter set in in the centre, slow dimerisation on the left and fast 
operator binding on the right. In each case, we divide the collection of TSE configurations 
according to the state of the operator. For each operator state, we plot histograms for 
the A-coordinate, in such a way that the area under a histogram corresponds to the total 
number of TSE points with that operator state. The histograms are colour coded according 
to operator state. Considering first the central panel of the upper row (baseline parameter 
set), we see that the green histogram (OA2) is shifted towards larger values of A than the 
red histogram (OB2) — i.e. the state of the operator and A are correlated in the TSE. 
This means that if a -B dimer is bound to the operator, then, on average, the number of A 
molecules has to exceed the number of B molecules in order to have the same value of Pb, 
and vice versa. We also see that the area under the OB2 histogram is larger than that under 
the OA2 histogram — indicating that the TSE has predominantly B2 bound to the operator, 
even though the switch is symmetric. Turning next to the right panel — rapid operator 
association and dissociation — we see that again the OA2 histogram is shifted towards larger 
values of A relative to the OB2. However, in this case, the areas under the two histograms are 
approximately equal. Thus increasing the rate of operator binding appears to have caused 
the transition state for switch flipping to become symmetric in A and B. The left panel shows 
the results for slow dimerisation, kf = 0.1. This plot is virtually indistinguishable from the 
baseline parameter results — indicating that changing the dimerisation rate has little effect 
on the transition state ensemble, as suggested by Fig. [6]A. These results unambiguously 
demonstrate that, for the exclusive switch, fluctuations arising from TF-DNA association- 
dissociation reactions are central to the flipping mechanism, while those arising from TF-TF 
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association-dissociation reactions have little effect on the ffipping mechanism, although they 
can inffuence the dynamics of the flipping trajectories and hence the switching rate. 

Drawing together the observations of Figs. [2], |4]A., [5]A., [6]A. and [71 we can now understand 
the dependence of the exclusive switch flipping rate on the rate of operator binding (Fig. 
[2B). In the limit of slow operator binding and unbinding 0, 0], the binding of the minor- 
ity species to the operator strongly enhances the flipping of the switch: when the minority 
species happens to bind the operator, it will stay on the DNA for a relatively long time, thus 
blocking the synthesis of the majority species and allowing the production of the minority 
species. In this limit, the system can reach the dividing surface with relatively few produc- 
tion/degradation events. As the rate of operator binding and unbinding is increased, each 
transition involves many operator binding/unbinding events, and consequently proteins of 
both species are produced and decay during the transition process. Here, the state of the 
operator is increasingly slaved to the difference in the total number of A and B molecules, 
A. In the adiabatic limit of fast operator binding, the probability that a molecule of type A 
or B is bound to the operator is completely determined by A In this limit, the dividing 
surface is located at A ~ and {Nq^^) ^ {Nqb^)', to reach the separatrix, the system has to 
wait for a series of fluctuations in the birth and decay of both species that lead to nA ~ wb- 
This implies that the total number of copies of A and B at the dividing surface decreases 
as the rate of operator binding increases (Fig. [5]A). Because a series of production/decay 
events are required to reach the separatrix, the probability p{q*) is decreased as the rate of 
operator binding increases (Fig. [2p). In addition, having reached the separatrix, the system 
requires more production/decay events to take it to the B state. This increases the proba- 
bility that it will "recross" the separatrix and eventually return to A instead of contributing 
to B — resulting in a decrease in the kinetic prefactor R in Fig. [2F. 

Figs. [5H7] suggest that the rate of dimerisation only has a marginal effect on the switching 
pathways. However, our view of the switching pathways naturally depends on the repre- 
sentation in which we choose to plot them. We have investigated many representations to 
see whether the rate of dimerisation could affect the switching pathways. Perhaps the most 
important one is the average number of dimers {NB2) as a function of {Nb{Nb — 1)). How- 
ever, also in this representation the rate of dimerisation only has a very minor effect on the 
switching pathways; in fact, near the top of the dividing surface, the dimerisation reaction 
is in steady state (data not shown). This supports our confusion that dimerisation affects 
the rate at which the transition paths traverse state space (and hence R), but not the route 
they take (and thus not p{q*)). 

The effect of TF-TF association/dissociation fluctuations on the dynamics of the trajec- 
tories can perhaps be understood by considering that in order to start a switching event 
from one stable state to the other, two copies of the minority species must be produced. 
They must then dimerise and bind to the operator, to shut down production of the ma- 
jority species. If the dimerisation rate is comparable to the degradation rate, it becomes 
increasingly probable that copies of the minority species, once produced, are removed from 
the system before they can form a dimer. Thus, decreasing the dimerisation rate actually 
reduces the chance that the switch can flip. This effect is truly dynamical in origin. We 
note that it is also fundamentally different from enhanced switch stability via cooperativity 
due to nonlinear degradation [28 . 

Lastly, while operator binding is an equilibrium reaction, it couples to reactions that 
are out of equilibrium, such as protein production and decay. As a result, the dynamics of 
operator binding can lead the exclusive switch to behavior that is unique for non-equilibrium 
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systems. This can be seen by comparing the forward paths from A to B with the backward 
paths from B to A in Fig. [6]A.. When the rate of operator binding is fast, the forward and 
backward paths essentially coincide. This situation differs markedly for the system with the 
base-line parameter set and for the system with slow dimerisation: although the switch is 
symmetric on interchanging A and B, the transition path ensemble (TPE) for the transition 
from A to B does not coincide with that from B to A j2|. This is a manifestation of the 
fact that this switch is a non-equilibrium system: for equilibrium systems that obey detailed 
balance and microscopic reversibility, the forward and backward paths must necessarily 
coincide. The fact that the forward and backward paths do not coincide also means that 
the switching paths do not follow the path of highest steady-state phase space density, 
which, for equilibrium systems, would correspond to the lowest free-energy path: Since this 
system is symmetric, this "lowest-free energy path" is symmetric on interchanging species 
A and B, while Fig. [6JA. shows that the dynamical switching trajectories are not (unless 
operator binding is fast). This also means that for this system it is essential not to make 
the Markovian assumption of memory loss, which underlies path sampling schemes such as 
Milestoning [29] and PPTIS 



C. Switching pathways - General switch 

We now turn our attention to the switching pathways for the general switch, again ob- 
tained with FFS. Fig. 03 shows for the three different parameter sets the switching tra- 
jectories as averaged in the Pb ensemble and projected onto the nA,nB plane. As for the 
exclusive switch, the forward and backward paths do not coincide, which, as mentioned 
above, reflects the fact that the genetic switch is a non-equilibrium system. However, in 
many other respects the pathways of the general switch differ markedly from those of the 
exclusive switch. Firstly, the switching trajectories of the general switch cross the dividing 
surface A = at very low values of ua and n-Q — on average, there is only one dimer of each 
species at the transition surface. Moreover, the paths display a sharp deviation when they 
reach the dividing surface. Lastly, paths obtained for different values of the rate constants 
essentially coincide (the black, red and green lines overlap). This last observation suggests 
that the transition paths are rather insensitive to variations in the rate constants of dimeri- 
sation and operator binding, an observation that should be contrasted with the observation 
that both p{q*) and R do depend upon the magnitude of those rate constants (see Fig. [3]). 

Fig. [6]B shows the state of the operator for every value of A during the transition, for the 
baseline parameter set (the curves for the other parameter sets are virtually indistinghuis- 
able). Initially, when the system is still in the basin of attraction of the stable state A, the 
operator is mostly in state OA2. However, as the system leaves this basin, the state of the 
operator rapidly becomes dominated by OA2B2. Indeed, this operator state, which is absent 
in the exclusive switch, plays a pivotal role in flipping the general switch. Its occupancy 
peaks at A ~ —5, corresponding to the top of the barrier that separates the stable state A at 
A = —27 from the metastable state at A = 0. Here, at A = 0, the occupation statistics of the 
operator is given by the equilibrium distribution [OA2]:[OB2]:[OA2B2], with [OA2] = [OB2]. 
As Fig. [33 shows, the transition state ensemble coincides with the metastable state around 
A = 0, and in this ensemble the operator is predominantly in state [OA2B2]. As the sys- 
tem leaves the metastable state towards the B state, the state of the operator progressively 
moves toward [OB2]. 

We are now able to explain the process of flipping the general switch. When a dimer of 
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the minority species is produced, it immediately binds to the operator and drives it in the 
state OA2B2. In this state, the production of both proteins is suppressed, and the system 
is depleted of almost all its components [13, The approach to the transition state is 



then driven mostly by a decrease of the majority species via protein decay rather than an 
increase of the minority species via protein production, as in the exclusive switch; this is 
the reason why the general switch crosses the diving surface at lower values of ua and hb 
than the exclusive switch. Importantly, if one of the two dimers leaves the operator, it can 
immediately rebind, thereby restoring the previous situation and allowing the transition to 
continue. By contrast, if the minority species leaves the operator in the exclusive switch, then 
most likely the majority species will bind the operator, thereby blocking further progress of 
the transition. This explains why both the pathways and the rate of flipping are much more 
sensitive to the rate of operator binding in the exclusive switch than in the general switch. 

The presence of the state OA2B2 also underlies the metastability of the general switch 
at A = (Fig. HP). As long as both species are present in the system, the state OA2B2 
is the most stable operator state, and in this state no proteins can be produced. As a 
consequence, a small fluctuation in A away from A = via the unbinding of, say, dimer A 
leading to the production of protein B, is not sufficient to flip the switch: most likely the 
dimer will rebind the operator, blocking further production of B; only when the dimer A 
dissociates and one of its monomers is degraded will the system commit itself to the stable 
state B. The probability that the dimer is degraded before it rebinds the operator increases 
as the rate of dimer dissociation increases; this is the origin of the increase of k^B, R and 
p{q*) with increasing dimersiation rate fcf for low fcf seen in Fig. El Finally, we note that 
the discrete character of the components in combination with their low copy number is 
important 31, 3^: flipping the switch away from the metastable state at A = requires 
the unbinding and subsequent degradation of essentially one molecule. The metastability 
is indeed absent in a mean-fleld continuum analysis that ignores the discrete nature of the 
components. 



V. DISCUSSION 

In this paper, we have analysed the stability and switch flipping dynamics of two types of 
bistable genetic toggle switches, as a function of the rates of transcription factor dimerisation 
and operator binding. This allows us to assess the influence of two key sources of fluctuations 
in the network on the overall system behaviour. 

We have varied the rate constants of the TF-TF and TF-DNA association/dissociation 
reactions over more than three orders of magnitude (see Figs. [2] and [3]) • This reflects the 
wide range of observed rate constants for cellular biochemical reactions. For instance, in 
prokaryotic cells, the inverse rate of protein production, k~^^^, is in the range seconds to 
minutes jssj. Since the size of a typical prokaryote is about Ipw? (based on E. coU), this 
corresponds to kp^odV = 10~^ — 10nM~^/min. The rate of monomer-monomer association, 
fcf, is about 10~^ — lO^^nM^^/min, while the dimer dissociation rate is of the order of 
kd = 10~^ — 10^/min, corresponding to dissociation constants in the range = — lO^nM 



28( 1 . This means that kf = 10^^ — lO/cprodV^- Figs. [2JA. and[3]A show that the switching rate 
kAB is fairly insensitive to changes in the dimerisation rate when kf > k^j-o^V, but is highly 
sensitive to dimerisation rate for kf < kp^o^V. This shows that the rate of dimerisation can 
strongly affect the network stability under biologically relevant conditions. Rate constants 
for protein-DNA association/dissociation are observed to vary over a similarly broad range 



13 



33l | . Figs. [2|3 and [3]B demonstrate that this variation can have a marked effect on ffipping 
rates for multistable networks in hving cells. 

The steady state phase space density in the region of the stable states is very robust 
to every parameter change (Fig. H]). Yet, changing the rate constants does strongly affect 
the switching between these states. Factorising the switching rate into the probability p{q*) 
of finding the system at the dividing surface, and a kinetic prefactor R, we find different 
results for the two versions of the switch: while for the exclusive switch dimerisation affects 
the switching rate predominantly via the kinetic prefactor, for the general switch it affects 
both the kinetic prefactor and the probability of being at the separatrix; on the other hand, 
operator binding affects the flipping rate of the exclusive switch both via R and p{q*), 
whereas the effects on the flipping rate of the general switch are exerted predominantly 
through a modification of its steady state distribution near the separatrix. 

These results can be understood by analyzing the transition paths and the transition state 
ensemble (TSE). This shows that, in the case of the exclusive switch, changes in the operator 
binding rate strongly affect the properties of the TSE, while the dimerisation rate has little 
effect on the TSE. We conclude that, for the exclusive switch, operator binding fluctuations 
play a crucial part in the reaction coordinate, while dimerisation fluctuations can affect the 
dynamics of the transition but have little effect on the route that it takes in phase space. The 
case of the general switch is rather different: here, the presence of the sterile, doubly-bound 
state OA2B2 makes the flipping pathways rather insensitive to both sources of fluctuations, 
even though the latter do affect the flipping rate. The resolution of this paradox lies in the 
fact that the switch is a non-equilibrium system: in contrast to equilibrium systems that 
obey detailed balance and microscopic reversibility, in non-equilibrium systems the forward 
and backward transition paths can form a cycle, as observed here; changing microscopic 
transition rates (i.e., reaction rate constants) can then change the stationary distribution 
p{q) and the flipping rate, even though the location of the transition paths in state space 
is unaltered. Protein-protein and protein-DNA association and dissociation reactions are a 
common feature of a wide range of biological control networks. We therefore hope that our 
results will be useful to understand the factors governing stability in multistable biochemical 
networks in general. 

Genetic switch flipping is an example of a non-equilibrium rare event. Rather few studies 
have been made of rare events in non-equilibrium systems, but a variety of simulation and 
analytical approaches have been developed to analyse rare events in equilibrium systems. 
Here, it is often assumed that one coordinate, the "reaction coordinate", is slow, while 
the other degrees of freedom are fast. In this case, the transition can be modelled by 
assuming that the reaction coordinate evolves according to a Langevin equation, while the 
other degrees of freedom play the role of friction. Although the concept of free energy is 
not applicable to non-equilibrium systems, one can nevertheless define a "barrier" which 
corresponds to the maximum of — logp(g), as we do in this paper. The results presented 
here show that that "barrier crossing" in the model toggle switch differs fundamentally 
from this classical scenario. For the genetic switch, the reaction coordinate consists of 
at least two parameters, namely the difference in total copy number of species A and B 
and the state of the operator [3]. Moreover, these coordinates evolve on comparable time 
scales — the operator state fluctuates on time scales similar to those of protein production 
and decay; in addition, their dynamics mix in a non-equilibrium fashion |9| — the degradation 
and production of proteins are non-equilibrium processes. This hampers the application of 
standard theoretical tools to model barrier crossings |9|]- New theoretical approaches may 
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be required to model such rare events in non-equilibrium systems. 
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LIST OF FIGURES 

1. Pictorial representation of the model switch, corresponding to reaction scheme ([2]). 
Two divergently-transcribed genes are under the control of a shared regulatory binding 
site on the DNA (the operator site O). Proteins A and B can bind, in homodimer 
form, to the operator. Each TF acts to block the production of the other species. In 
the exclusive switch, only one type of TF can bind at any given time (meaning that 
never the production of both species can be suppressed), whereas in the general switch 
both types of TF can bind (in which case the production of both species is repressed). 

2. Panels A and B show the switching rate kAB for the exclusive switch as a function of 
the dimerisation rate k{ (A) and the rate of operator binding kon (B). Dissociation rates 
are scaled such that the equilibrium constants remain constant: = k\,/ki = 1/V 
and = k^s/kon = 1/{5V). Panels C and D show the probability p{q*) of being at 
the dividing surface, as a function of kf (C) and kon (D). Panels E and F show the 
kinetic prefactor, as defined by Eq. ([1]), as a function of k{ (E) and fcon (F). 

3. Panels A and B show the switching rate k^B for the general switch as a function of the 
dimerisation rate kf (A) and the rate of operator binding fcon (B). Dissociation rates 
are scaled such that the equilibrium constants remain constant: = k^^/kf = 1/V 
and = kos/kon = 1/(51^). Panels C and D show the probability p(g*) of being at 
the dividing surface, as a function of kf (C) and kon (D). Panels E and F show the 
kinetic prefactor, as defined by Eq. ([1]), as a function of /cf (E) and kon (F). 

4. Probability distribution as a function of the order parameter A = UA — riB, with nx the 
total copy number of species X, for the exclusive switch (A) and for the general switch 
(B). The distributions are obtained with FFS calculations [lH, for three different sets 
of parameters. 

5. Switching paths projected onto the uaiUb surface, for three different sets of parame- 
ters. (A) Paths averaged in the Pb ensemble for the exclusive switch, where ua and 
ub are averaged over configurations with the same value of Pb, where nx is the total 
copy number of species X. The forward paths, corresponding to transitions from A to 
B are shown with solid lines, while the reverse transitions, from B io A are shown 
with dashed lines. (B) Switching paths for the general switch. In this projection, the 
paths are highly insensitive to variations in parameters. The hypersurface A = is 
crossed with a lower total number of A and B molecules. 

6. A) Exclusive switch: probability that a B2 dimer is bound to the operator, {Nq-q^), as a 
function of Pb for three different sets of parameters. The solid lines correspond to the 
transition from A to 5, while the dashed lines corresponds to the reverse transition 
from B to A. B) General switch: operator occupancies during the transition from 
A io B and vice versa (the empty state O is not shown since it is always scarcely 
occupied), for the baseline parameter set; the results for the other parameter sets are 
indistinghuisable. 

7. The probability p{\) for the transition state ensemble {Pb = 0.5) for the transition 
from A to B^ where A = ua — ub- Top row A) Exclusive switch; Bottom row B) 
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General switch. The probabihty p(A) is spht into colour-coded contributions from 
the three operator states; the area under each histogram gives the probability (Nqx) 
that the operator is bound to species X (the three areas thus sum to unity). The 
left panels correspond to the system with slow dimerisation fcf = 0.1; the middle 
panels correspond to the system with the base-line parameters; the panels on the right 
corresponds to the system with fast operator binding kon — 500. 
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FIG. 1: Morelli et al. 
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